真菌基因组分析报告

真菌基因组分析报告

1 项目概况

1.1 项目信息

合同编号 BN-200515NJ01S02N3
合同名称 江苏省中国科学院植物研究所一个真菌全基因组图谱开发
甲方单位 中国科学院植物研究所
样品名称 Cg201
联系邮箱
撰写时间 2022年01月11日

1.2 项目背景

基因组类型：真菌基因组

预计基因组大小：40 M-60 M

1.3 项目完成情况

组装得到的基因组长为62,777,491 bp，其中包含28 条contig，平均contig长度为2,242,053.25 bp，最长的contig为7,950,879 bp，GC含量为50.12 %。

2 生物信息学分析

对测序得到的序列进行低质量数据过滤得到有效数据，基于过滤后数据，进行组装及后续分析，最终得到完整的基因组序列。生物信息分析流程见下图：

图 1 生物信息分析流程图

注：绿色部分为个性化分析，不包含在标准流程内。

2.1 测序数据及其质量控制

2.1.1 测序数据格式

高通量测序得到的原始图像数据文件经碱基识别 （Base Calling）分析转化为原始数据（raw reads 或 rawdata），结果以 FASTQ[1] （简称为 fq）文件格式存储，其中包含测序序列（reads）的碱基信息以及其对应的测序质量信息。Raw Data 以 FASTQ 格式存储，每个 read 由四行描述，如下：

表 1 FASTQ 格式示例

@ST-E00600:415:HFV5LALXX:1:1101:21486:1031 1:N:0:TGACAGGT
GATCTGGTCCTCTCGGACCGCGATTCTTGCCTTCTGTGAACCTCGACAGATCGATCTTGAACCTCGACAGATCGATCTTGAACCTCGA
+
AAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAFFFAAF

注：第一行以 “@” 开头，为 Illumina 测序标识符（Sequence Identifie）和描述文字（选择性部分），该行相当于 reads ID；第二行是碱基序列；第三行以 “+” 开头，后为 Illumina 测序标识符，一般为节省存储空间，省略 “+” 后测序标识符；第四行是对应碱基的测序质量，该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。

2.1.2 测序质量值

在测序过程中会存在一定的错误率，测序质量值[2]可衡量碱基未正确检出的概率。如果测序错误率用 P 表示，碱基质量值用 Q_score 表示，则有如下关系：

\[Q_{-score} = -10 \times log_{10}P\]

每个碱基测序错误率是通过测序质量值（Q_score） 通过公式转化得到，而测序错误率是在碱基识别（Base Calling）过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所示：

表 2 测序错误率与 Phred score 对应表

Phred Quality Score Probability of Incorrect Base Call Base Call Accuracy
Q10 1/10 90%
Q20 1/100 99%
Q30 1/1000 99.9%
Q40 1/10000 99.99%

注：测序错误率分布具有两个特点：（1）测序错误率会随着测序序列（Sequenced Reads）长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致错误；（2）前几个碱基的位置会发生较高的测序错误率，推测前几个碱基测序错误率较高的原因是由于随机引物碱基导致。

碱基质量值越高表明碱基识别越可靠，碱基测序错误的可能性越小。如，对于碱基质量值为 Q20 的碱基识别，100 个碱基中有1个会识别出错；对于碱基质量值为 Q30 的碱基识别，1,000 个碱基中有 1 个会识别出错；Q40 表示 10,000 个碱基中才有 1 个会识别出错。

2.1.3 数据过滤标准

由于原始测序数据可能包含低质量序列、接头序列等，为了保证信息分析结果的可靠性，原始测序数据需要经过过滤，从而得到 clean reads，仍然以 FASTQ 格式存储（该序列文件可直接用于发表，公共数据库提交等，后续分析都基于此数据）。

二代过滤标准：

1） 去除N碱基含量超过 5% 的 reads； 2） 去除低质量（质量值小于等于 5）碱基数目达到 50% 的 reads； 3） 去除有 adapter 污染的 reads； 4） 去除 PCR 扩增造成的重复序列。

三代过滤标准：

1） 去除平均质量值小于等于 7 的序列。

2.1.4 测序数据量统计

2.1.4.1 三代测序数据量统计

测序得到的所有数据中，去除 fail reads（平均质量值小于等于 7 的序列）后，得到 pass reads（即平均质量值大于7的序列），用于后续分析。

表 3 三代测序数据量信息统计

Rank Flag TotalBase TotalReads MaxLen AvgLen N50 L50 N90 L90 meanQ
>0 all 8,439,038,676 547,220 178,810 15,421.65 29,108 102,233 7,848 297,875 8.61
>0 pass 7,547,769,304 482,378 178,810 15,647 29,242 91,136 7,981 264,539 8.92
>0 fail 891,269,372 64,842 171,130 13,745.24 28,002 11,115 6,769 33,479 6.29
>5000 all 7,932,448,328 351,642 178,810 22,558.3 30,467 93,727 11,626 250,663 8.66
>5000 pass 7,105,823,496 313,961 178,810 22,632.82 30,569 83,746 11,682 223,762 8.94
>5000 fail 826,624,832 37,681 171,130 21,937.44 29,580 9,991 11,205 26,913 6.32
>10000 all 7,342,680,716 269,494 178,810 27,246.17 32,038 84,290 15,465 211,329 8.66
>10000 pass 6,582,332,237 241,037 178,810 27,308.38 32,120 75,395 15,504 188,972 8.94
>10000 fail 760,348,479 28,457 171,130 26,719.2 31,283 8,902 15,154 22,360 6.33
>50000 all 1,076,320,367 18,322 178,810 58,744.69 57,390 8,219 51,159 16,194 8.70
>50000 pass 972,154,475 16,558 178,810 58,712.07 57,384 7,434 51,171 14,636 8.94
>50000 fail 104,165,892 1,764 171,130 59,050.95 57,501 785 51,055 1,558 6.38
>100000 all 9,585,398 86 178,810 111,458.11 108,423 40 100,832 77 8.22
>100000 pass 7,577,634 68 178,810 111,435.79 109,252 32 100,806 61 8.69
>100000 fail 2,007,764 18 171,130 111,542.44 105,898 9 102,634 17 6.47

注：Rank 为数据长度梯度，>0 即为所有数据；Flag 为数据类型，all 为全部测序数据，pass 为有效测序数据，fail 为过滤掉的数据；TotalBase 为碱基数目；TotalReads 为 reads 条数；MaxLen 为数据的最大长度；AvgLen 为数据的平均长度；N50 为数据的 N50，将所有的 reads 按照从长到短排序后依次累加，当相加的长度达到所有 reads 总长度的一半时，最后一个加上的 read 长度即为 N50；L50 为数据的 L50，将所有 reads 按照从长到短排序后依次累加，当相加的长度达到所有 reads 总长度的一半的总序列条数；N90 为 reads 的 N90，算法同 N50；L90 为数据的 L90，算法同 L50；meanQ 为平均质量值。

图 2 测序数据的 reads 长度和 reads 条数分布图

注：横坐标为测序的 reads 长度，纵坐标为不同长度对应的 reads 条数，两条黑竖线从左到右分别为所有数据的平均长度和 N50，不同颜色分别对应 pass reads 和 fail reads。

图 3 测序数据的 reads 长度和数据量分布图

注：横坐标为测序的 reads 长度，纵坐标为不同长度对应的数据量，两条黑竖线从左到右分别为所有数据的平均长度和 N50，不同颜色分别对应 pass reads 和 fail reads。

图 4 测序数据的质量值和 reads 条数分布图

注：横坐标为每个 read 的平均质量值，纵坐标为不同测序质量值对应的 reads 数目，不同颜色分别对应 pass reads 和 fail reads。

2.1.4.2 二代测序数据量统计

二代测序，原始测序数据为 rawdata，经过滤后，得到有效测序数据（cleandata）。

表 4 二代测序数据量信息统计

Sample_name Total_reads Total_bases GC_content Q20 Q30
rawdata 32,831,468 4,924,720,200 52.07% 97.93% 94.24%
cleandata 32,831,236 4,924,685,400 52.07% 97.93% 94.24%

注：Type 为数据类型；Total_reads 为测序 reads 数； Total_bases 为测序总碱基数；GC_Content 为 G/C 碱基数占总碱基数量的百分比；Q20、Q30 分别为 Phred 数值大于 20、30 的碱基占总体碱基的百分比。

2.2 K-mer 分析

利用测序得到的 reads，采用基于 K-mer 的分析方法来估计基因组大小和杂合率等情况。K-mer 是指一段长度为 K bp 的序列。从一段连续序列中迭代地选取长度为 K 个碱基的序列，若序列的长度为 L，K-mer 长度为 K，那么可以得到 L-K+1 个 K-mer。我们对测序得到的 reads 取 K-mer，然后统计每个 K-mer 出现的频数。根据 Lander_waterman 算法，基因组大小（G）满足如下公式：

\[C_{base} = C_{k-mer} \times \frac{L}{L-K+1}\]

\[G =\frac{n_{k-mer}}{C_{k-mer}}=\frac{n_{base}}{C_{base}}\]

注：Cbase 和 Ck-mer 为覆盖碱基的期望深度和 K-mer 期望覆盖深度，nbase 和 nk-mer 为序列的碱基总数和 K-mer 总数。在数据量一定的情况下，K-mer 的深度频数是服从泊松分布的，因此 K-mer 深度频数分布的峰值为对应的深度，作为 K-mer 期望深度的估计值。

K-mer 分布图大致分为以下几种：（1）当只有一个主峰时，该个体为纯合体或单倍体；（2）在 x = a 处出现主峰，x = 2a 处有一次峰，说明一部分片段出现的期望值是大部分的 2 倍，这些片段为重复片段，次峰为重复峰；（3）在 x = a 处出现主峰，x = 0.5a 处有一个次峰，说明部分片段出现的期望值是大部分的 1/2，当序列有杂合时，包含杂合位点的 K-mer 因为分成了两部分，所以出现频率变为一半，次峰为杂合峰；（4）当出现两个主峰，且峰高相差不大，两峰横坐标又是 2 倍关系，说明该个体高杂合或高重复。

图 5 Kmer-Depth 和 KmerSpecies - Frequence 分布图

注：蓝色线表示实际的 K-mer 曲线，黑线是用模型估计的 k-mer 曲线，黄线是 unique 数据对应的 K-mer 曲线，红线表示由于测序错误导致的错误曲线，虚线表示推测的 K-mer 峰值。

使用genomescope（version 1.0.0）软件对基因组大小进行估计，基因组的评估结果见下表：

表 5 基因组大小估计

Item Number
Heterozygosity 0.064%
Genome Haploid Length 59,486,073 bp
Genome Repeat Length 13,321,670 bp
Genome Unique Length 46,164,404 bp
Model Fit 98.406%
Read Error Rate 0.019%

注：Heterozygosity 为杂合度，Genome Haploid Length 为基因组大小，Genome Repeat Length 为重复序列的长度，Genome Unique Length 为非重复序列的长度，Model Fit 为模型匹配度，Read Error Rate 为 read 错误率。

2.3 基因组组装

使用 NECAT 软件进行基因组纠错和拼接，得到初始拼接结果；然后使用 Racon（version：1.4.11）软件基于三代测序数据对拼接结果进行两轮纠错后，再进行两轮二代 reads 的 Pilon（version：1.23）纠错，最后使用 purge_haplotigs 对纠错后的基因组进行去杂合后，得到最终的组装结果。

表 6 组装结果数据统计

Item Value
Total_length(bp) 62,777,491
Total_length_withoutN(bp) 62,777,491
Total_number 28
GC_content(%) 50.12
N50(bp) 5,448,682
N90(bp) 4,249,034
Average(bp) 2,242,053.25
Median(bp) 437,801.00
Min(bp) 36,391
Max(bp) 7,950,879

注：Total_length 为组装长度；Total_length_withoutN 为组装结果中不含 gap 的长度；Total_number 为组装的序列条数；GC_content 为 GC 含量；N50 为数据的 N50，将所有的序列按照从长到短排序后依次累加，当相加的长度达到总长度的一半时，最后一个加上的序列长度即为 N50；N90 算法同 N50；Average 为平均长度；Median 为中位长度；Min 为最小长度；Max 为最大长度。

2.3.1 基因组组装评价

基因组组装评价一般从以下几个方面：二代比对率、覆盖度以及 BUSCO 评价。

2.3.1.1 二代数据比对

二代比对率和覆盖度一般大于 90%，认为基因组组装较好。

表 7 二代数据比对统计

Type Number
map_rate 99.62%
Average_depth 77.95
Coverage 97.55%

注：map_rate 为二代数据比对率；Average_depth 为平均覆盖深度；Coverage 为覆盖率。

2.3.1.2 BUSCO 评价

使用 BUSCO 软件（version：4.1.4）基于真菌数据库（ fungi_odb10 ）来评估基因组组装的完整性，BUSCO 评价一般大于 90%，认为基因组组装较好。

表 8 BUSCO 评价统计

Item Number Percent(%)
Complete BUSCOs (C) 288 99.3
Complete and single-copy BUSCOs (S) 287 99.0
Complete and duplicated BUSCOs (D) 1 0.3
Fragmented BUSCOs (F) 1 0.3
Missing BUSCOs (M) 1 0.3
Total BUSCO groups searched 290 100.0

注: Complete BUSCOs 为完整的 BUSCO；Complete and single-copy BUSCOs 为完整且单拷贝的 BUSCO；Complete and duplicated BUSCOs 为完整但多拷贝的 BUSCO；Fragmented BUSCOs 为不完整的 BUSCO；Missing BUSCOs 为缺失的 BUSCO；Total BUSCO groups searched 为参考的 BUSCO 数目。

2.4 基因组结构预测

通过基因结构预测，能够获得基因组详细的基因分布和结构信息，也将为功能注释和进化分析工作提供重要的原料。

2.4.1 编码基因预测

编码基因预测包括预测基因组中的翻译起始位点和终止位点、内含子和外显子区域以及蛋白质编码序列等。基因预测主要使用 BRAKER 软件（version：2.1.4），首先使用 GeneMark-EX 训练模型后，调用 AUGUSTUS 进行预测。

表 9 编码基因注释统计

Type Number
the total number of gene 15,845
the average of mRNA_length 1,759.90
the average of cds_length 1,492.18
the average of exon_number 3.08
the average of exon_length 483.71
the average of intron_length 128.07
the total number of exon 58,684
the total number of intron 39,661
the total intron length 5,079,443

注：the total number of gene 为总基因数；the average of mRNA_length 为平均的 mRNA 长度；the average of cds_length 为平均的 CDS 长度；the average of exon_number 为每基因包含的平均外显子数；the average of exon_length 为平均的外显子长度；the average of intron_length 为平均的内含子长度；the total number of exon 为总外显子数；the total number of intron 为总内含子数；the total intron length 为总的内含子长度。

2.4.1.1 编码基因的 BUSCO 评价

对预测的基因用 BUSCO 软件（version: 4.1.4）基于真菌数据库（ fungi_odb10 ）进行完整性评估。统计结果图如下：

图 6 编码基因 BUSCO 评价结果图

注：Single-copy 为单拷贝的 BUSCO；duplicated 为多拷贝的 BUSCO；Fragmented 为片段化的 BUSCO；Missing 为缺失的 BUSCO。

2.4.2 非编码 RNA 注释

非编码 RNA，指的是不被翻译成蛋白质的 RNA，这些 RNA 都具有重要的生物学功能。其中 tRNA 、rRNA 直接参与蛋白质的合成。使用 INFERNAL（version：1.1.2）基于 Rfam 数据库进行各类 ncRNA 预测，并分类统计。

表 10 非编码 RNA 注释结果统计

Class number totalLen(bp) meanLen(bp)
rRNA 215 415,690 1,933
sRNA 3 751 250
snRNA 32 4,379 136
tRNA 377 31,481 83

注：rRNA 为核糖体 RNA；tRNA 为转运 RNA；sRNA 为小的调控 RNA；snRNA 为核仁小 RNA。

2.4.3 重复序列注释

重复序列可分为散在重复序列和串联重复序列。散在重复序列又称转座子元件，包括4种：LTR、LINE、SINE 和 DNA 转座子。根据重复多少可以分为高度重复序列、中度重复序列和低度重复序列。使用 RepeatModeler 软件（version：1.0.4）构建自身的 repeat 库，合并 repbase 库后，使用 RepeatMasker（version：4.0.5）进行基因组的重复序列注释。

表 11 重复序列注释结果统计

Item Subfamily Number Length(bp) Coverage
SINE / 14 942 0.00%
LINE / 683 62,306 0.10%
LTR / 1,929 374,610 0.60%
LTR Gypsy 1,000 250,921 0.40%
LTR Copia 733 110,577 0.18%
DNA / 932 81,990 0.13%
Satellite / 139 15,146 0.02%
Simple_repeat / 17,494 792,180 1.26%
Low_complexity / 1,649 82,606 0.13%
Other / 66 5,131 0.01%
Unknown / 39 9,184 0.01%
Total / 22,945 1,404,379 2.24%

注：SINE 为短散在元件；LINE 为长散在元件；LTR 为长末端重复，主要包括两种类型，Gypsy 和 Copia；DNA 为转座子；Satellite 为卫星重复序列；Simple_repeat 为简单重复；Low_complexity 为低复杂性重复；Other 为其他类型重复；Unknown 为未知的重复序列；Total 为总的重复序列。

2.4.4 次级代谢基因簇分析

次级代谢产物是微生物在一定的生长时期，以初级代谢产物为前体合成的对微生物的生命活动无明确功能，并非生长繁殖所必需的物质。采用antiSMASH程序（version6.0.0）对基因组进行预测。

表 12 次级代谢基因簇预测统计表

2.4.4.1 Overview页面注释

报告对Overview和cluster做了简要说明，更多详细说明见官方文档，官方链接为：https://docs.antismash.secondarymetabolites.org/understanding_output/

注：第一列：cluster名称，对应于结果文件夹中region0gbk；

第二列：antiSMASH自定义的一些缩写来表示不同类型的次级代谢产物簇；

第三、四列：cluster的起始终止位置；

第五列：最相似的结果；

第六列：已知代谢簇；

第七列：相似性值。

2.4.4.2 cluster页面注释

注：Overview所在行：1.1表示基因组第1个染色体上的第1个cluster，1.2表示基因组第1个染色体上的第2个cluster；以此类推；条形图一二行显示基因簇的生物合成类型和位置；第三行条形图是基因是按预测函数编码的颜色。假定生物合成基因为红色，运输相关基因为蓝色，调控相关基因为绿色。

示例图注

2.5 基因功能注释

2.5.1 通用数据库注释

基因的功能注释是指根据已有数据库对基因功能和所参与的代谢途径进行标注，包括 Motif、结构域、蛋白质功能及所在的代谢通路等信息的预测。为获得全面的基因功能信息，进行了九大数据库的基因功能注释，包括：Nr、Pfam、eggCOG、Uniprot、KEGG、GO、Pathway、Refseq、Interproscan。

基因功能注释主要有如下 2 种方法：（1）序列相似性搜索：将全长基因中基因编码的蛋白序列与现有蛋白质数据库 Uniprot、Refseq、NR以及代谢通路数据库 KEGG 进行 diamond blastp（version ：2.9.0；参数：–evalue 1e-5）比对，获得序列的功能信息，以及蛋白可能参与的代谢通路信息。其中 KEGG 注释使用 KOBAS（version: 3.0）关联到 KEGG ORTHOLOGY 以及 PATHWAY。Uniprot 数据库记录了每个蛋白质家族与 Gene Ontology 中的功能节点的对应关系，通过此系统预测基因编码的蛋白序列所执行的生物学功能。基于数据库（Uniprot/Swiss-Port）之间的关联，得到 eggNOG 数据库的注释信息，并挑选 COG 注释结果，进行 COG 的分类统计及绘图。

（2）Motif相似性搜索：蛋白质中，一般由一个或多个功能区构成，这些区通常被称为域。结构域的不同组合方式产生了多种多样的蛋白质。因此蛋白结构域的鉴别对分析蛋白质的功能来说尤其重要。使用 hmmscan（版本：3.1；参数：e-value 0.01）进行结构域预测，获得蛋白质的保守序列、motif和结构域等。Pfam 数据库是一个蛋白质家族大集合，依赖于多序列比对和隐马尔可夫模型。利用 InterProScan（version：5.33-72.0）对二级数据库 InterPro 子数据中的 CDD-3.16、Coils-2.2.1、Gene3D-4.2.0、Hamap-2018_03、MobiDBLite-2.0、Pfam-32.0、PIRSF-3.02、PRINTS-42.0、ProDom-2006.1、ProSitePatterns-2018_02、ProSiteProfiles-2018_02、SFLD-4、SMART-7.1、SUPERFAMILY-1.75以及 TIGRFAM-15.0 数据库进行比对，获得蛋白质的保守序列、motif 和结构域等。

表 13 注释统计表

Item Count Percentage
All 19,016 100%
Annotation 19,016 100%
Uniprot 10,038 52.79%
Pfam 14,427 75.87%
Refseq 8,555 44.99%
Nr 18,838 99.06%
Interproscan 14,154 74.43%
GO 9,916 52.15%
KEGG 5,272 27.72%
Pathway 3,284 17.27%
COG 1,592 8.37%

注：Annotation 为至少有一种注释的基因；Uniprot 为注释到 Uniprot 数据库的基因；Pfam 为注释到 Pfam 数据库的基因；Refseq 为 Refseq 数据库注释到的基因；Nr 为注释到 Nr 数据库的基因；Interproscan 为注释到 Interproscan 数据库的基因；GO 为注释到 GO 数据库的基因；KEGG 为注释到 KEGG 数据库的基因； Pathway 为注释到 KEGG Pathway 数据库的基因；COG 为注释到 COG 数据库的基因。

2.5.1.1 基因功能查询

对 9 大数据库的注释信息进行合并，方便对基因的搜索和挖掘。

表 14 注释汇总表

注：gene_id 为基因 ID；Uniprot 为 Uniprot 注释；Pfam 为 Pfam 注释；Refseq 为 Refseq 注释；Nr 为 Nr 注释； Interproscan 为 Interproscan 注释；COG 为 COG 注释；KEGG 为 KEGG 注释； Pathway 为 Pathway 注释；GO 为 GO 注释。

2.5.1.2 NR 注释结果统计

通过 Nr 库比对注释的结果，统计并绘图，结果见下图。根据比对结果，统计比对最多的前 10 个物种，其余划分到 Other species 类，物种分布图如下图所示。

图 7 Nr 数据库比对到序列的物种分布图

注：Nr 注释到的不同的同源物种分布情况。不同颜色代表不同的物种。

2.5.1.3 Gene Ontology 分类

将 GO（Gene Ontology）注释信息经简化后得到 GOslim 分类，将基因的功能从细胞组分、分子功能、生物学过程三个方面进行汇总统计后，选取每个分类下前 20 个注释最多的 GOslim 的二级分类进行绘图，结果如下图示。

图 8 GO 功能注释分类统计图

注：横坐标为 GO 各分类内容，纵坐标为基因数目。此图展示的是在全部基因背景下 GO 各二级功能的基因富集的情况，体现此背景下各二级功能的地位。

2.5.1.4 COG 分类

COG（Cluster of Orthologous Groups of proteins）数据库是基于细菌、藻类、真核生物的系统进化关系构建得到的数据库，利用 COG 数据库可以对基因进行直系同源分类。COG 分为 26个 group, 将 COG 注释的基因按 COG 的 group 进行分类，结果如下图所示。

图 9 COG 功能注释分类统计图

注：横坐标为 COG 各分类内容，纵坐标为基因数目。在不同的功能类中，基因所占多少反映对应时期和环境下代谢或者生理偏向等内容，可以结合研究对象在各个功能类别的分布作出科学的解释。

2.5.1.5 Pfam 结构域分类

基于 Pfam 结构域的注释，对每个结构域注释的基因进行统计汇总，并将注释最多的前 20 个结构域绘图展示，结果如下图所示。

图 10 Pfam 功能注释分类统计图

注：Pfam 数据库是一系列蛋白质家族的集合，其中每一个蛋白家族都以多序列比对和隐马尔科夫模型的形式来表示。图中横坐标是蛋白质家族的名称，纵坐标是比对到该蛋白质家族的基因的数。

2.5.1.6 KEGG Pathway 分类

对序列做 KEGG 注释后，根据它们参与的 KEGG 代谢通路进行分类，结果如下图所示。

图 11 KEGG Pathway 功能分类图

注：横坐标为 Pathway 分类下注释的基因数目；纵坐标为 Pathway 分类，不同颜色代表所属的不同的大分类。

2.5.2 专有数据库注释

2.5.2.1 碳水化合物酶注释

碳水化合物在很多生物学功能中具有重要地位，通过研究碳水化合物相关酶可以得到大量有意义的生物学信息，CAZy 数据专注于分析碳水化合物酶的基因组 、结构和生物化学信息。用 HMMER（version：3.2.1，过滤参数 E-value < 1e-18；coverage > 0.35）基于 CAZy 数据库对蛋白质序列进行注释。

表 15 碳水化合物活性酶数据库功能分类对应关系表

功能 ID
糖苷水解酶 （Glycoside Hydrolases ，GHs）
糖基转移酶 （Glycosyl Transferases，GTs）
多糖裂合酶 （Polysaccharide Lyases，PLs）
碳水化合物酯酶 （Carbohydrate Esterases，CEs）
辅助氧化还原酶 (Auxiliary Activities , AAs)
碳水化合物结合模块 （Carbohydrate-Binding Modules，CBMs）

表 16 CAZy 注释

Family HMM HMM length Query ID Query length E-value(how similar to the family HMM) HMM start HMM end Query start Query end Coverage
NA 548 g10003.t1 548 0 8 547 26 546 0.9835766
NA 307 g10016.t1 363 0 3 273 46 301 0.8794788
NA 342 g10029.t1 920 0 1 342 592 917 0.9970760
NA 239 g10045.t1 393 0 2 238 51 291 0.9874477
NA 693 g10077.t1 748 0 9 599 21 669 0.8513709
NA 130 g10098.t1 263 0 4 129 64 195 0.9615385

注：Family HMM 为蛋白质序列名称；HMM length 为模型蛋白的长度；Query ID 为待检测的基因 ID；Query length 为待检测的基因长度；E-value 为比对的 E-value 值，该值越低越好。

2.5.2.2 病原宿主互作分析

PHI 全称为 Pathogen Host Interactions，即病原与宿主互作数据库。该数据库收录的内容经过实验验证，主要来源于真菌、卵菌和细菌病原感染的宿主包括动 物、植物、真菌以及昆虫。使用 Diamond blastp （version >：2.9.0；参数：–evalue 1e-5）基于 PHI 数据库对目的蛋白质序列进行注释。

表 17 PHI 注释

qseqid sseqid pident length mismatch gapopen
g5.t1 I1RXV6 37.468 395 188 12
g5.t1 I1RXV6 26.957 230 142 4
g5.t1 Q6DQW3 44.484 281 141 8
g5.t1 Q6DQW3 29.369 412 245 11
g5.t1 Q6DQW3 25.495 404 237 13
g5.t1 Q0UI00 32.292 384 222 11

注：qseqid 为检测序列；sseqid 为比对上数据库的目标序列；pident 为比对的相似度；length 为比对长度；mismatch 为比对的错配碱基数目；gapopen 为比对开的 gap 数目；(结果文件剩余列项含义： qstart 为检测序列比对的起始位置；qend 为检测序列比对的终止位置；sstart 为比对上数据库序列的起始位置；send 为比对上数据库序列的终止位置；evalue 为比对的评价指标；bitscore 为比对的得分；PHI_info 为PHI的注释信息）。

2.5.2.3 耐药基因预测

通过 CARD 数据库注释，获得每个基因组中包含的耐药基因的信息。 CARD 数据库是以 Antibiotic Resistance Ontology (ARO) 为分类单位的形式所构建，用于关联抗生素模块及其目标、抗性机制、基因变异等信息。比对结果可 展示每一个在 CARD 数据库中注释到的基因对应的位置，ARO 编号以及 ARO 分类描述等，通过 ARO 分类描述可了解每个基因与抗生素抗性相关的具体功能。

表 18 CARD 数据库注释

ORF_ID Best_Hit_ARO ARO AMR Gene Family
g1.t1 IMI-4 3001861 IMI beta-lactamase
g2.t1 vanL 3002910 glycopeptide resistance gene cluster; van ligase
g2.t2 vanL 3002910 glycopeptide resistance gene cluster; van ligase
g3.t1 mecC 3001209 methicillin resistant PBP2
g5.t1 lmrB 3002813 ATP-binding cassette (ABC) antibiotic efflux pump
g6.t1 iri 3002884 rifampin monooxygenase
g7.t1 lsaA 3000300 ABC-F ATP-binding cassette ribosomal protection protein
g8.t1 vanWI 3003724 vanW; glycopeptide resistance gene cluster
g9.t1 mef(B) 3003107 major facilitator superfamily (MFS) antibiotic efflux pump
g10.t1 tet(41) 3000569 major facilitator superfamily (MFS) antibiotic efflux pump

注：ORF_ID 为待检测基因 ID；Best_Hit_ARO 为比对最好的 ARO；ARO 为耐药基因的 ID；AMR Gene Family 为耐药基因家族。

2.5.2.4 细胞色素 P450 注释

细胞色素 P450（Cytochromes P450，简称 CYP450）是一大类以亚铁血红素为辅基的蛋白家族。它们能催化许多种底物的氧化反应。因这类蛋白的还原态与 CO 结合后 会在 450 nm 处检测到最大吸收波长，故又命名为 P450。它参与内源性物质和包括药物、环境化合物在内的外源性物质的代谢。使用 Diamond blastp （version >：2.9.0；参数：–evalue 1e-5）基于 FungalP450 数据库对目的蛋白质序列进行注释。

表 19 细胞色素 P450 注释统计表

Database Number
CYP450 2,141

2.5.2.5 毒力基因预测

DFVF 数据库，全称 database of fungal virulence factors，是专门用于研究真菌毒力因子的数据库。利用Diamond blastp （version ：2.9.0；参数：–evalue 1e-5）将预测得到的蛋白序列与毒力因子数据库（DFVF）进行比对，分析测序菌株中存在的毒力相关基因。

表 20 毒力基因预测

Gene_ID VFDB_ID Description
g10011.t1 SUB1_ARTOC Subtilisin-like protease 1
g10480.t1 LAP2_ASPFU Probable leucine aminopeptidase 2
g10600.t1 CUTI1_COLGL Cutinase 1
g10666.t1 SUB1_TRIRU Subtilisin-like protease 1
g10893.t2 BOT2_BOTFU Presilphiperfolan-8-beta-ol synthase
g10893.t3 BOT2_BOTFU Presilphiperfolan-8-beta-ol synthase

注：Gene_ID 为基因 ID；VFDB_ID 为毒力因子 ID; Description 为毒力因子为功能描述。

2.5.2.6 TCDB 注释

TCDB 是对膜转运蛋白（Membrane Transport Protein）进行分类的一个数据库，它制定了一套转运蛋白分类系统（Transporter Classification）, 简称 TC System, 类似于对酶进行分类的 EC 系统，只不过 TC 系统同时提供了功能和进化信息；TCDB 对于每一个转运蛋白家族，提供了一个 TC Nmuber, TC Number 由小数点分隔的5为数字或者字母构成。

表 21 TCDB 注释统计表

Database Number
TCDB 2,820

2.5.2.7 信号肽预测

使用软件 SignalP（version:5.0）对所有的预测到的基因的蛋白序列进行分析，找出含有信号肽的蛋白，预测结果如下：

表 22 信号肽预测统计表

Protein_type Number
Signal peptide protein 2,443
图 12 信号肽预测示例图

注：横坐标为蛋白序列；纵坐标为信号肽出现的可能性。

2.5.2.8 跨膜蛋白预测

使用软件 tmhmm（version:5.0）所有的预测到的基因的蛋白序列进行分析，找出含有跨膜螺旋的蛋白，即为跨膜蛋白，预测结果如下：

表 23 跨膜蛋白预测统计表

Protein_type Number
Transmembrane protein 4,211

2.5.2.9 分泌蛋白预测

在上述预测的含有信号肽的蛋白中去除含有跨膜螺旋的蛋白，剩余的蛋白即为分泌蛋白，预测结果如下：

表 24 分泌蛋白预测统计表

Protein_type Number
Secreted protein 2,444

3 实验流程

3.1 样品检测

1） 样品的外观是否含有异物； 2） 琼脂糖电泳检测是否样品有降解以及 DNA 片段大小； 3） Nanodrop 检测 DNA 纯度； 4） Qubit 精确地对 DNA 进行定量。

3.2 二代实验流程

3.2.1 二代测序文库构建及库检

DNA 样品检测合格后，使用 Covaris 超声波破碎仪随机打断，再经末端修复、加 A 尾、加测序接头、纯化、PCR 扩增等步骤完成整个文库制备工作。文库构建完成后，先使用 Qubit 2.0 进行初步定量，稀释文>库，随后使用 Agilent 2100 对文库的插入片段进行检测，插入片段大小符合预期后，使用 Q-PCR 方法对文库的有效浓度进行准确定量，以保证文库质量。

图 13 二代建库流程

3.2.2 Illumina 上机测序

文库检测合格后，按照有效浓度及目标下机数据量的需求将不同文库 pooling 至 Flow cell，cBOT 成簇后使用 Illumina 高通量测序平台Illumina NovaSeq 进行测序。

3.3 三代实验流程

3.3.1 三代测序文库构建以及质量检测

样本质检合格后，用 Megaruptor 对基因组 DNA 进行随机打断；利用磁珠富集、纯化大片段 DNA；使用 BluePippin 全自动核酸回收仪将大片段进行切胶回收；对片段化的 DNA 进行损伤修复；纯化后在 DNA 片>段两端进行末端修复并进行加 A，使用 SQK-LSK109 连接试剂盒中接头进行连接反应，最后用 Qubit 精确地对建好的 DNA 文库进行定量检测。

图 14 三代文库构建流程

3.3.2 Nanopore 上机测序

建库完成后将一定浓度和体积的 DNA 文库加入到 Flow cell 中，并将 Flow cell 转移到 Oxford Nanopore PromethION 测序仪进行实时单分子测序。

4 软件及数据库列表

表 25 软件及数据库汇总表

软件/数据库 版本 分析 网址
genomescope 1.0.0 基因组survey https://github.com/schatzlab/genomescope
NECAT 20200119 三代组装 https://github.com/xiaochuanle/NECAT
nextDenovo 2.3.1 三代组装 https://github.com/Nextomics/NextDenovo
Racon 1.4.13 三代序列纠错 https://github.com/isovic/racon
Purge Haplotigs Latest Version 去除杂合重复 https://bitbucket.org/mroachawri/purge_haplotigs
minimap2 2.17-r941 长序列比对 https://github.com/lh3/minimap2
bwa 0.7.17 短序列比对 https://github.com/lh3/bwa
samtools 1.10 比对结果处理 https://github.com/samtools/samtools
pilon 1.23 二代序列纠错 https://github.com/broadinstitute/pilon
purge_haplotigs Latest Version 去除杂合重复 https://github.com/skingan/purge_haplotigs_multiBAM
braker 2.1.5 编码基因注释 https://github.com/Gaius-Augustus/BRAKER
BUSCO 4.1.4 完整性评估 https://github.com/RoyNexus/busco
INFERNAL 1.1.3 ncRNA注释 https://github.com/EddyRivasLab/infernal
RepeatModeler 1.0.4 构建重复序列模型 https://github.com/Dfam-consortium/RepeatModeler
RepeatMasker 4.0.7 重复序列鉴定与分类 https://github.com/rmhubley/RepeatMasker
BLAST+ 2.9.0 比对 https://github.com/enormandeau/ncbi_blast_tutorial
Diamond 2.0.4 比对 https://github.com/bbuchfink/diamond
KOBAS 3.0 kegg注释 https://github.com/xmao/kobas
HMMER 3.2.1 pfam注释 https://github.com/EddyRivasLab/hmmer
SignalP 5.0 信号肽预测 http://www.cbs.dtu.dk/services/SignalP/
tmhmm 2.0 跨膜蛋白预测 http://www.cbs.dtu.dk/services/TMHMM/
InterProScan 5.33-72.0 InterProScan注释 https://github.com/ebi-pf-team/interproscan
Uniport 2019_07 Uniport功能注释 https://www.uniprot.org/
refseq 89 refseq功能注释 https://www.ncbi.nlm.nih.gov/refseq/
Pfam 31.0 结构域注释 https://pfam.xfam.org/
GO 2014 GO功能注释 http://geneontology.org/
KEGG 87.0-20180701 KEGG功能注释 https://www.kegg.jp/kegg/
COG 2014 直系同源基因注释 https://www.ncbi.nlm.nih.gov/COG/
Nr 201909 Nr注释 https://ftp.ncbi.nlm.nih.gov
Rfam Latest Version ncRNA注释 http://rfam.xfam.org/
CAZy Latest Version 碳水化合物活性酶分析 http://www.cazy.org/
PHI Latest Version 病原菌与宿主互作分析 http://www.phi-base.org
CARD Latest Version 耐药基因注释 https://card.mcmaster.ca/
CYP450 Latest Version 细胞色素P450注释 http://drnelson.utmem.edu/CytochromeP450.html
DFVF Latest Version 毒力基因预测 http://sysbio.unl.edu/DFVF/
TCDB Latest Version 膜转运蛋白预测 http://www.tcdb.org/
antiSMASH 6.0.0 次级代谢基因簇分析 https://docs.antismash.secondarymetabolites.org/understanding_output/

5 联系我们

客 服 电 话：027-62435310
客 服 邮 箱：service@benagen.com
公 司 网 站：www.benagen.com
公 司 地 址：武汉东湖新技术开发区高新大道888号高农生物园总部B区12C栋

6 参考文献

[1] Cock P J, Fields C J, Goto N, et al. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants[J]. Nucleic Acids Research, 2010, 38(6):1767.
[2] Ewing B, Green P. Base-Calling of Automated Sequencer Traces Using Phred. II. Error Probabilities[J]. Genome Research, 1998, 8(3):186-94.
[3] Liu B, Shi Y, Yuan J, et al. Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects[J]. Quantitative Biology, 2013, 35(s 1–3):62-67.Ruan, J. and Li, H. (2019) Fast and accurate long-read assembly with wtdbg2. bioRxiv. doi:10.1101/530972
[4] Walker B J, Abeel T, Shea T, et al. Pilon: An Integrated Tool for Comprehensive Microbial Variant Detection and Genome Assembly Improvement[J]. PLOS ONE, 2014, 9.
[5] Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform[M]. 2009.
[6] Li H E A. The Sequence Alignment / Map (SAM) Format[J]. Bioinformatics, 2009, 25(1 Pt 2):1653-4.
[7] Haas B J, Delcher A L, Mount S M, et al. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies[J]. Nucleic Acids Research, 2003, 31(19):5654-5666.
[8] Stanke M, Diekhans M, Baertsch R, et al. Using native and syntenically mapped cDNA alignments to improve de novo gene finding[J]. Bioinformatics, 2008, 24(5):637-644.
[9] Simao F A, Waterhouse R M, Ioannidis P, et al. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[J]. Bioinformatics, 2015: btv351.
[10] Katharina J. Hoff, Simone Lange, Alexandre Lomsadze, Mark Borodovsky and Mario Stanke, “BRAKER1: Unsupervised RNA-Seq-Based Genome Annotation with GeneMark-ET and AUGUSTUS”, Bioinformatics, 2015, Nov 11
[11] Ter-Hovhannisyan et al Gene prediction in novel fungal genomes using an ab initio algorithm with unsupervised training, Genome Research (2008) 18, pp 1979-1090
[12] Mario Stanke, Mark Diekhans, Robert Baertsch, David Haussler (2008). Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics, 24(5), pages 637–644, doi: 10.1093/bioinformatics/btn013
[13] Büttner, E., Gebauer, A. M., Hofrichter, M., Liers, C., & Kellner, H. (2019). Draft Genome Sequence of Xylaria longipes DSM 107183, a Saprotrophic Ascomycete Colonizing Hardwood. Microbiology resource announcements, 8(12), e00157-19. doi:10.1128/MRA.00157-19
[14] G. Gremme, V. Brendel, M.E. Sparks, and S. Kurtz. Engineering a software tool for gene structure prediction in higher organisms. Information and Software Technology, 47(15):965-978, 2005
[15] Benjamin Buchfink, Chao Xie & Daniel H. Huson (2015). Fast and Sensitive Protein Alignment using DIAMOND. Nature Methods, 12, 59–60 (2015) doi:10.1038/nmeth.3176. PMID:25402007
[16] Ai, C., Kong, L. CGPS: A machine learning-based approach integrating multiple gene set analysis tools for better prioritization of biologically relevant pathways.
[17] Jones P, Binns D, Chang H Y, et al. InterProScan 5: genome-scale protein function classification[J]. Bioinformatics, 2014, 30(9):1236-1240.
[18] Simon C Potter, Aurélien Luciani, Sean R Eddy, Youngmi Park, Rodrigo Lopez, Robert D Finn, HMMER web server: 2018 update, Nucleic Acids Research, Volume 46, Issue W1, 2 July 2018, Pages W200–W204, https:/ /doi.org/10.1093/nar/gky448
[19] Lowe, T.M. and Chan, P.P. (2016) tRNAscan-SE On-line: Search and Contextual Analysis of Transfer RNA Genes. Nucl. Acids Res. 44:W54-57.
[20] Lagesen K, Hallin PF, Rødland E, Stærfeldt HH, Rognes T Ussery DW RNammer: consistent annotation of rRNA genes in genomic sequences, Nucleic Acids Res. 2007 Apr 22.